Applying Adaptive Time-Dependent DMRG to Calculate 
the Conductance of Strongly Correlated Nanostructures 
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A procedure based on the recently developed "adaptive" time-dependent density-matrix-renormalization- 
group (DMRG) technique is presented to calculate the zero temperature conductance of nanostructures, such 
as a quantum dots (QD's) or molecular conductors, when represented by a small number of active levels. The 
leads are modeled using non-interacting tight-binding Hamiltonians. The ground state at time zero is calcu- 
lated at zero bias. Then, a small bias is applied between the two leads, the wave-function is DMRG evolved in 
time, and currents are measured as a function of time. Typically, the current is expected to present periodici- 
ties over long times, involving intermediate well-defined plateaus that resemble steady states. The conductance 
can be obtained from those steady-state-like currents. To test this approach, several cases of interacting and 
non-interacting systems have been studied. Our results show excellent agreement with exact results in the non- 
interacting case. More importantly, the technique also reproduces quantitatively well-established results for 
the conductance and local density-of-states in both the cases of one and two coupled interacting QD's. The 
technique also works at finite bias voltages, and it can be extended to include interactions in the leads. 

PACS numbers: 73.63.-b,71.27.+a,72.10-d,85.65.+h 
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I. INTRODUCTION 

The rapidly developing investigations in the area of 
nanometer-scale systems and its concomitant potential tech- 
nological applications in real devices have induced consider- 
able interest in the study of electrical transport through small 
molecules and quantum dots. In fact, the construction of 
molecular electronic devices** is among the most exciting 
areas of investigations in physics, and theoretical guidance 
is needed for the success of this vast effort. Molecules can 
change their shape and position relative to the leads as elec- 
trons enter or leave the molecule, making the study of these 
systems very challenging. Moreover, Coulomb correlations 
cannot be neglected in small devices. For a conceptual under- 
standing of these complex systems, it is imperative to develop 
models and unbiased many-body methods that rely on a mini- 
mal number of assumptions, in order to accurately handle both 
strong Coulombic and electron-phonon couplings. Quantum 
dots constructed using conventional semiconductor technol- 
ogy also provide an important playground for the analysis of 
transport properties in nanoscopic systems, and the theoretical 
challenges in this context are equally important 3 . 

The conductance of small nanoscopic systems can be the- 
oretically estimated using a variety of techniques. Among 
the most popular approaches are the ab-initio calculations us- 
ing density functional theory (DFT). These one-electron self- 
consistent methods have been successful in describing various 
I-V characteristics 4 *** However, the applicability of these 
ideas must be carefully scrutinized, as recently remarked by 
Muralidharan et al? . For example, it is clear that in small 
molecules charging effects are important, and they effectively 
act as quantum dots in the Coulomb Blockade regime. More- 
over, techniques that do not take into account the strong corre- 
lation between electrons cannot capture important effects such 



as the Kondo resonance (arising from the coupling between 
localized spins and conduction electrons) which induces a 
new channel for transport in a variety of small systems**. In 
addition, it is well known that several bulk materials, such 
as transition metal oxides, cannot be described with ab-initio 
methods that neglect correlations. The complexity of their be- 
havior, including potentially useful effects such as large mag- 
netoresistances in Mn oxides**, may manifest in nanoscopic 
systems as well, and the use of strongly correlated materials 
in nanodevices may lead to interesting applications. To study 
all these systems (small molecules, quantum dots, and in gen- 
eral nanodevices that include strongly correlated materials), 
techniques beyond DFT must be developed. 

A similar challenge occurred before in the study of bulk 
materials and several years of research have shown that the 
following two-steps process leads to profound insights. The 
first step consists of a simple modeling of the material, typ- 
ically either deducing the relevant degrees of freedom from 
atomistic considerations when the states are very localized, or 
borrowing from band structure calculations to isolate the min- 
imal ingredients needed to capture the essence of the problem. 
The second step, the hardest, is solving the resulting model, 
which is typically of a tight-binding nature with the addition 
of Coulombic and phononic couplings. In the strong coupling 
regime, the use of numerical techniques provides the most re- 
liable unbiased approach for the approximate investigation of 
tight-binding-like models with Coulombic interactions**. As 
a consequence, a natural path toward the study of transport 
in strongly correlated nanoscopic systems can also start with 
simple models and use computational techniques for their 
analysis. It is the main purpose of this paper to propose a tech- 
nique that can be used to study transport in systems described 
by strongly correlated electronic models. 

The method to calculate conductances proposed here re- 
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lies on the successful Density Matrix Renormalization Group 
(DMRG) technique 1213 . While in principle this method is 
applicable only to quasi-one-dimensional models, such a ge- 
ometry is quite acceptable in several nanoscopic important 
problems including transport in arrays of quantum dots or us- 
ing small molecules as bridges between leads. However, a 
straightforward application of the original DMRG methodol- 
ogy is not immediately useful to study transport, particularly 
when arbitrarily large external electric and magnetic fields, 
with potentially complicated time dependences, are switched 
on and off at particular times. Nevertheless, progress toward a 
computational tool for these type of problems has been steady 
in recent years. For instance, important numerical methods 
for dynamical DMRG studies were presented to handle fre- 
quency dependent spectral functions 14 . More directly focus- 
ing on real-time investigations, interesting techniques were 
proposed 1 ^. While useful for many qualitative applications, 
these methods are in general not as accurate and stable as 
needed for the detailed study of finite bias transport in compli- 
cated nanosystems. The reason is that the method of Re f Ha is 
"static" in the sense that the truncated Hilbert space found to 
be optimal at time t=Q, namely before switching on the exter- 
nal fields, is kept at all times. This approach breaks down after 
relatively short times, since extra states are needed for a proper 
description of transport at finite times. A further approxima- 
tion to improve on this first proposal is to enlarge substantially 
the initial Hilbert space so that it remains suitable for proper- 
ties calculated at finite timesA£>. This technique has the prob- 
lem that the number of states grows rapidly with the simula- 
tion time and eventually it becomes impractical. Nevertheless, 
the method has been successfully used to study the propaga- 
tion of a density excitation in an interacting clean systemii. 

Recently, important developments have led to the "adap- 
tive" time-dependent version of the DMRG method, that is 
efficient over long times and, thus, it is suitable to handle 
the problems we are focusing on (for a detailed review see 
Ref 18). The method to be used here was developed indepen- 
dently by White and Feiguir. 19 and Daley et al*&, after the 
idea of how to do time-evolution to a matrix product was in- 
troduced by Vida! 21 , and relies on an adaptive optimal Hilbert 
space that follows the state as time grows. The method is 
based on a Suzuki-Trotter break-up of the evolution opera- 
tor, and as a consequence a Trotter truncation error is intro- 
duced. Fortunately, this systematic error can be easily esti- 
mated and controlled. The adaptive DMRG numerical method 
will only be briefly reviewed below since our proposal uses 
the technique to calculate conductances, but does not mod- 
ify the method itself. The reader should consult the origi- 
nal literatur e) 19 i 20 for more details. It is important to remark 
that the technique is easy to implement once a ground-state 
DMRG code is prepared and, moreover, the time-evolution is 
stable, as shown explicitly in our results and in some previous 
investigations (further improvements can be added with the 
time-step targeting method recently proposed by Feiguin and 
White 22 ). The conclusion of our effort documented below is 
that the adaptive method provides accurate results for the cal- 
culation of conductances. The technique has passed the test 
of noninteracting electrons, as well as the cases of one and 



two interacting quantum dots, where a subtle Kondo effect 
occurs. Moreover, the method is not restricted to small bi- 
ases but it produces reasonable answers at finite bias as well. 
As a consequence, it has the potential of being the method of 
choice to study transport under both weak and strong external 
fields, in small nanostructures of substantial complexity. Mul- 
tilevel model Hamiltonians, possibly inspired by ab-initio cal- 
culations, can be used to describe the "bridge" between leads. 
In addition, the method is particularly transparent since it re- 
lies on the straightforward calculation of a current in the pres- 
ence of a voltage, rather than relying on other indirect linear- 
response formulas. 

Of course, the reader must be aware that the method is not 
of unlimited applicability. If the molecules or Kondo clouds 
are too long in size, eventually not even DMRG can handle the 
very long chains needed for a proper description. For com- 
pleteness, and to assure a balanced description of the tech- 
nique, one of these difficult cases is also presented in our 
manuscript. But often the qualitative physics can be under- 
stood by relaxing parameters, thus we expect that even in very 
complicated cases the propose technique will be helpful, at 
least at the conceptual level. Other limitations of the present 
technique is that energy dissipation is not incorporated, and 
the temperature is restricted to be zero. Improving on these 
issues is a task left for the near future. 

It is important to remark that there are other numerical tech- 
niques that can also be used to study transport in strongly cor- 
related nano-systems. One of them is the Numerical Renor- 
malization Group, that evolved from Wilson's original RG 
ideas. This technique is quite accurate, as exemplified by 
some recent calculations 23 , but it cannot be used for arbitrary 
problems. Since our goal is to try to develop a method that 
can handle the fairly complex models that will be used in the 
near future to represent, e.g., small molecular conductors, this 
method does not have sufficient flexibility for our purposes. In 
cases where NRG works, it should be the method of choice, 
but this occurs in a small subset of problems in the area of 
transport across correlated systems. A second approach relies 
on the static DMRG method, using a ring geometry and with 
a current induced by a flux threading the ring2i. A recently 
proposed third method combines linear response Kubo theory 
with static DMRG and the conductance is calculated based on 
correlation functions in the ground stata 2 ^. It wo uld be in- 
teresting to find out if the methods of Refs[24 25 can reach 
a similar accuracy as ours for the case of the one and two 
quantum dots. A fourth method is the Exact Diagonalization 
technique followed by a Dyson equation embedding proce- 
dure (ED+DE) 26,27 where the interacting region is solved ex- 
actly (including some sites of the leads), and then the rest of 
the leads are taken into account via a Dyson self-consistent 
approach. The method directly treats bulk systems, contrary 
to the DMRG technique that is necessarily limited to a large 
but finite chain, it is flexible and has led to interesting results 
for difficult cases, such as center-of-mass phonons in molec- 
ular conductors, and multilevel systems 28,29,30,31-32 . However, 
the Dyson embedding is somewhat arbitrary and it is difficult 
to control its accuracy. It is our intention in the near future 
to combine the ED+DE method with the DMRG approach 
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the left (right) lead. £f c iustcr is the Hamiltonian of the clus- 
ter where the interactions are present. Finally, Hciuster-ieads 
is the Hamiltonian that connects the interacting region to the 
leads, typically a hopping term. 



FIG. 1 : Schematic representation of the geometry used in our study. 
The leads are modeled by tight-binding Hamiltonians. The ground 
state at time zero is calculated at zero bias. Then, a finite bias AV is 
applied between the two leads (without ramping time for simplicity, 
but this can be changed in future studies) and the resulting current is 
measured. 



discussed in this paper, and for the physical problems where 
these independent techniques give similar results, then the 
case can be made that a reliable conclusion has been reached. 
Thus, ED+DE and the present method are complementary. 

The organization of the manuscript is as follows. After the 
present introduction, in Section II, the models are defined and 
the technique is very briefly described. Section III contains the 
important test of noninteracting electrons (note that although 
the Coulombic coupling is zero, there are different hopping 
amplitudes at different links). Here, the systematic behavior 
of the method is discussed in detail. Section IV deals with 
the case of a quantum dot, with a nonzero Hubbard coupling. 
The value of U is comparable to the hoppings, to prevent the 
Kondo cloud from reaching huge sizes that would render the 
DMRG method useless (nevertheless, one "large" U case is 
studied for completeness, to illustrate the limitations of the 
method). Section V contains the case of two dots, which 
themselves can be coupled in a "singlet" preventing conduc- 
tion or loosely coupled having individually a Kondo effect. 
The conclusions in Section VI briefly summarize our findings. 

II. MODEL AND CONDUCTANCE CALCULATION 

In general, the systems studied here consist of a relatively 
small region where Coulomb interactions are present, weakly 
coupled to two non-interacting leads (see Fig^. The interact- 
ing region can represent one or several quantum dots (QD's), 
a single-molecule conductor, or other nanoscopic regions. In 
fact, the generality of the method presented in this paper, al- 
lows for a wide variety of interacting systems. 

The leads are modeled as ideal tight-binding chains. As ex- 
amples, the focus will be on one QD and two QD's connected 
in series. The total Hamiltonian of theses systems can be writ- 
ten in general as 

H = -ffloads + -^cluster + ^cluster— leads; (i) 

where H\ ca ds is the Hamiltonian of the leads, which is 

-Hleads = ~t leads /l c UtT C H+^ + c tia c n+lcr + k.C.]. (2) 

iieads 1S the hopping amplitude in the leads, which in the fol- 
lowing is taken as the energy scale (i.e. ii oa ds = 1)- The op- 
erator c'n (c ria ) creates an electron with spin a at site i in 



A. One Quantum Dot 

For the case of one QD, represented simply by one active 
level, -f/ciustcr can be written as 

"cluster VgUd + Undyridi, (3) 

where the first term represents the location of the energy level 
of the QD controlled by the gate voltage V g . The second term 
represents the Hubbard repulsion between electrons of oppo- 
site spins occupying the QD. rid = n-d^ + n di is the number 
of electrons at the dot. ff c luster— leads can be written as 

^c-lcads = -f'y^cjp.Crfg + c\ a Cda + fl.C.], (4) 
a 

where t' is the amplitude for the electronic hopping between 
the QD and the leads. ct CT creates and electron with spin a 
at the dot, while c] creates an electron at the last site of the 
left lead and c\. a creates an electron at the first site of the right 
lead, if sites are numbered from left to right. 

B. Two Coupled Quantum Dots in Series 

In the case of two QD's, //cluster can be written as 

^cluster = ^2 [ V g n da + Undain da i]-t''y^ j [c dl(T C d 2 t7 + h.C.], 
a=l,2 a 

(5) 

where rid a — ridot\ + n dai is the number of electrons at the 
quantum dot a, and t" is the hopping between the two dots. 

^cluster -leads is written as 

-^cluster -leads = ~t' ^Mo-Cfilff + 4<r c d2cr + h.C.}. (6) 



C. Conductance Calculation 

The current at any time t between nearest-neighbors sites i 
and j is calculated as 

Jij(t) = £)<*(t)l(4<*r - cj a Cv)|*(i)}, (7) 

a 

where is the wave function of the system at time t, 

which will be calculated with the DMRG method, using a 
number M of states in the process. c| CT creates an electron 
with spin a at site i, which can be part of the interacting re- 
gion or be at the leads. In the results presented below, the 
current shown without any link or site index corresponds to 

J(t) = (J L (t) + J R (t))/2, (8) 
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where Jl (t) is the current between the last site of the left lead 
and the first dot, and Jr{€) is the current between the last dot 
and the first site of the right lead, moving from left to right. 
The conductance G can be obtained by simply dividing the 
steady-state current by the total bias AV. The individual volt- 
ages ± A V/2 in the leads are applied uniformly in each one, as 
indicated in Fig[2 Note that the use of a symmetrized current 
J{t) is not crucial, since a good left-right symmetry is ob- 
served when currents at several links are studied, as described 
in more detail below. 



D. Technique 

Closely following Refill a brief description of the numer- 
ical technique is here provided for completeness. The basic 
idea is to incorporate the Suzuki-Trotter (ST) decomposition 
of the time-evolution operator 21 into the DMRG finite-system 
algorithm 19 - 20 . The second order ST decomposition of the ID 
Hamiltonian as employed in Refill can be written as 

e -irH ^ e -iTH 1 /2 e -iTH 2 /2^ e -irH 2 /2 e -iTH 1 /2^ ™ 

where Hj is the Hamiltonian of the link j. The DMRG repre- 
sentation of the wavefunction at a particular step j during the 
finite-system sweep is 

H') = ^ ]aj+1 r\l)\aj}\a j+ i)\r), (10) 

where I and r represent the states of the left and right blocks 
(in a truncated basis, optimally selected as eigenvectors of a 
density matrix), while on and oy+i represent the states of the 
two central sites. An operator A acting on sites j and j + 1 
(namely, only involving nearest-neighbors) can be applied to 
\ip) exactly, and re-expressed in the same optimal basis as 

[■A-1p]lctjOtj+ir = J ^-a j a j + 1 ;a'.a' j+1 ' l Pla' j a' j + 1 r- (H) 

a'. a'. . , 

Thus, the time evolution operator of the link j can be ap- 
plied exactly on the DMRG step j. As a consequence, the 
time evolution is done by applying e r lTB -^l 2 at DMRG step 
1, e ~ lTH2/2 at DMRG step 2, and so on, thus forming the 
usual left-to-right sweep. Then, applying all the reverse terms 
in the right-to-left sweep. A full sweep evolves the system one 
time step r. The error introduced by the second order decom- 
position is order r 3 in each time step 19 . Thus, upon evolving 
the system one time unit (1/t steps), an order t 2 error is in- 
troduced. Numerically, the influence of this small systematic 
error is easy to control. This brief summary gives the reader a 
rough idea of the technique used here. Details regarding lat- 
tice sizes, number of states kept in the DMRG procedure, and 
influence of other parameters are discussed below. 

III. NON-INTERACTING CASE 

Properties of the method discussed in this paper are exem- 
plified in Fig|2ja), where the current at the center of the chain 



is shown (divided by the voltage difference) as a function of 
time. The current in this figure is exactly calculated, not us- 
ing DMRG, since for non-interacting particles the problem 
reduces to a single electron problem. A small bias Ay=0.001 
will be used, unless otherwise stated. Thus, in this first study 
the focus will be on trying to reproduce results expected from 
linear response, but a few results with a finite bias are also 
included in the manuscript, as discussed below. Returning to 
Fig|2ja), for a bulk system the current would be expected to 
raise for a small fraction of time, and then reach a steady state. 
This indeed occurs even in our finite-size systems. In fact, the 
transition from zero current at t=0 to an approximately time- 
independent current regime is very fast, and can be barely ob- 
served in the scale of Fig|3Ja). But the existence of a very 
flat plateau in the current is clear, and its value will be used 
to extract the conductance below. Note that due to the finite 
size of the system, the current cannot continue in the same 
steady state at all times. The Hamiltonian is particle number 
conserving and, as a consequence, the presence of a current 
implies a population/depopulation of the leads, which cannot 
continue forever. In fact, once the front of the charge pulse 
reaches the end of the chain, it bounces back and eventually 
generates a current of the opposite sign. This effect will be 
discussed in more detail later. Here, it is important to remark 
that in spite of this periodicity present in finite open-end sys- 
tems, the flat plateaus are clearly defined over an extended 
period of time for the lattice of 402 sites used, and the value 
of the conductance can be easily deduced from those individ- 
ual plateaus, as discussed below. Note that the setup of Fig^ 
and the existence of plateaus in the current Fig|2]are natural 
in the DMRG/transport context and was observed before 33 . 
Our main contribution will be the use of the adaptive DMRG 
method for the calculations, as shown below. 

To further illustrate the propagation of charge in the cluster 
after the finite bias is switched on at time t=Q, in Fig|3]the 
exact current at different positions x is shown, parametric with 
time. At small times, t=5 (in units of tilt) only the central 
portion is affected as expected. At time £=55, the affected 
region is much larger, while at £=105, the front has reached 
the ends of the chain and soon after it starts bouncing back. 
At times £=200 and 205, approximately the initial condition is 
recovered, and almost everywhere the current to the left and 
right cancel out nearly exactly. For larger times, a reverse sign 
current is created. Note that in our studies there are no sources 
of dissipation, and the current will keep on oscillating forever. 
Adding inelastic processes is a next major challenge in this 
context, left for the future. 

It is also important to show that the existence of the plateaus 
is not restricted to very long chains of hundreds of sites, 
but they are visible on much smaller systems, increasing the 
chances that the numerical DMRG method can be used even 
for complicated nanosystems. Figure 0b) contains the cur- 
rent vs. time for a variety of lattice sizes, ranging from 402 to 
systems as small as 16 sites. The time width of the plateaus 
depends on L, as expected, but the value of the current at the 
plateaus is approximately L independent even up to systems 
as small as containing L=32 sites. Even the L=16 chain has 
a periodicity with a first plateau in the current which is also 
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FIG. 2: Exact results for J(t)/AV (in units of e 2 /h) vs. time (in 
units of ft/tioads), for the non-interacting 1QD case obtained with 
clusters of different lengths (L) and AV=0.001. (a) J(t)/AV ob- 
tained with a large cluster (L=402). J(t)/AV shows clear steady- 
state plateaus at ±2e 2 //i. The periodic changes in the current direc- 
tion are caused by its reflection at the open boundaries of the cluster, 
(b) J(t)/AV obtained with decreasing L. The steady-state plateau 
is obtained even with L = 32. The current is quasi-periodic with 
a period proportional to L. The parameters used are V s =U=0 and 
t'=0A. 



10 20 30 40 50 



10 



20 30 




/ — ^<L/ 






— Exact 






-- M = 200 






M= 100 






— M = 50 


(d)\ : 





15 20 



2 4 6 8 10 
t 



FIG. 4: DMRG results compared to the exact results for J(t)/AV 
obtained using different clusters L and number of states M, with 
Al/=0.001. (a) L = 96, (b) L = 64, (c) L = 32, and (d) L = 16. Note 
that for L = 96 and 64, M = 200 shows good qualitative agreement 
and M > 300 even shows good quantitative agreement with the exact 
results. For L = 32 and 16, M = 200 and 100 already show excellent 
quantitative agreement with the exact results. 
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FIG. 3: This figure shows the propagation of the current in the cluster 
after the bias AV=0.001 is applied at t=Q. The different panels show 
the current as a function of position x along the chain, at different 
times t. The size is L = 402 and the dot is at the site 201. The 
results were obtained exactly, since the Hubbard couplings are zero. 



in good agreement with the expected value from larger sizes. 
Thus, this behavior appears to be robust and the plateaus are 
also expected to be present for the chains that the DMRG 
method can handle. That this is the case can be shown in 
Fig|4] where DMRG results for the current vs. time are shown, 
compared with exact data. Consider first a sufficiently long 
chain, as shown in (a), such that a sharp plateau is observed 



in the exact result. Figure|^a) shows that increasing the num- 
ber of states M used in the DMRG approximation, a conver- 
gence to the exact solution is observed. In fact, for A/=300 
or higher, the DMRG results cannot be distinguished from the 
exact ones. A similar behavior is found using shorter chains as 
for the case with L=64 sites Fig|4jb), but in this example the 
plateau can be observed accurately even with a smaller num- 
ber of states such as M=200. The trend continues for smaller 
systems (c,d): For L=32 and 16, the DMRG method repro- 
duces the exact results with high accuracy using ~100 states. 

Figure |3 a) shows the conductance deduced from the be- 
havior of the current obtained with the DMRG method, versus 
the gate voltage V g , for the case of a single "noninteracting" 
quantum dot, namely one having U=Q. The hopping ampli- 
tude between the dot and the leads is i'=0.4. It is expected 
that the maximum value of the conductance be obtained when 
the level in the dot is aligned with the Fermi level of the leads, 
and this occurs in our case at T4=0. The DMRG results beau- 
tifully confirm this expectation. As the gate voltage changes 
away from 0, the conductance is expected to decrease sym- 
metrically and this is indeed shown in Fig[3a). In fact, the 
results at nonzero gate voltage are also in excellent quantita- 
tive agreement with the exact results. 

All the previous results were obtained for a sufficiently 
small value of the bias voltage A\7=0.001, as already ex- 
plained. It is interesting to observe how the results change 
when larger values of Ay are employed. Figure |5Jb) shows 
the current for a couple of biases. The existence of the plateau 
is clear in both cases, but for Ay=0.01, asymmetries between 
positive and negative bias can be observed, which are not ex- 
pected in the limit AV-^Q. As a consequence, the rest of 
the results discussed below were obtained with AV=0.001, 
unless otherwise noted. An easy criterion to realize if a suffi- 
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FIG. 5: (a) DMRG and exact results for G vs. V g in the case of 
one non-interacting QD, and t' = 0.4 (with AV=0.001). The DMRG 
results are obtained with L = 64 and M = 300. In this case, G is 
obtained from the value of the steady state current plateau. The exact 
results are for infinite leads. The results show a resonant tunneling 
peak of FWHM = At' 2 at V g = 0. (b) J(t) /AV for V g = ±0.3 and 
AV = 0.01, 0.001. Note that decreasing the bias voltage reduces 
the asymmetry between positive and negative V g and gives a better 
steady-state plateau. 



ciently small bias is used to obtain the linear response limit is 
to repeat the calculations for the same amplitude of bias, but 
opposite signs, and see if a noticeable difference is obtained. 

The method proposed here also works in the case of a finite 
bias voltage, namely it is not restricted to the linear response 
regime. To show that the technique can handle even a large 
bias, in Fig|5]results for the current vs. time are shown at the 
indicated voltages. It is only at AV as large as 1.0 that small 
differences are visually observed in the figure, between the 
DMRG and exact results. This can be easily fixed increasing 
slightly the number of states M, Thus, overall the method 
appears to be sufficiently robust to handle arbitrary voltages, 
showing the generality of the technique here proposed. Nev- 
ertheless, further work in the finite bias context will be im- 
portant to fully test this case, calculating differential conduc- 
tances and analyzing the regime of very strong bias. 

In our investigations, the numerical study was also carried 
out using a "static" procedure, where the t=0 DMRG basis 
is not expanded with growing time. In this case, the results 
are obtained by integrating the time-dependent Schrodinger 
equation using the fourth order Runge-Kutta method, and also 
using the DMRG ground state as the initia l state d This is 
to be contrasted with the procedure of Refs 19 20 where the 
basis is modified with time. Figure[7]shows the results of both 
procedures: clearly using an adaptive basis provides superior 
data, reproducing accurately the exact results. In the static 
procedure a similar accuracy is reached only by increasing 
substantially the number of states, thus missing its economical 
CPU-time advantages. 

The method discussed here also works nicely for the case 
of two non-interacting QD's, as shown in Fig|8]for two dif- 
ferent values of the hopping amplitude t" between the dots. 




FIG. 6: (a) DMRG results compared to the exact results in the case 
of one non-interacting QD for several intermediate and large values 
of AV, namely exploring the influence of a finite bias in the calcu- 
lations. The parameters used are Vg = U = 0, and t' = 0.4. Both 
DMRG and exact results are obtained with a cluster L = 32, using 
M = 200 states, for the DMRG results. 
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FIG. 7: The results obtained with the Suzuki-Trotter approach and 
the static Runge-Kutta method compared to the exact results for (a) 
AV = 0.001 and (b) AV = 0.5. The parameters used are V g = U = 
0.0, t' = 0.4, and L = 64. 



The slight difference between the DMRG and the bulk exact 
results can be improved increasing the number of sites. 

Although not directly related with the method to obtain 
the conductance of an interacting nano-system discussed in 
this paper, for completeness we have also studied the lo- 
cal density-of-states which are important to guide the intu- 
ition and contrast with other methods and scanning tunneling 
microscopy experiments as well. In Fig|9j results for non- 
interacting quantum dots are shown (namely, dots where the 
Hubbard repulsion is 0). Clearly, both the exact results (which 
are shown already in the bulk limit) and the DMRG results, 
slightly smeared by shifting, using a small imaginary compo- 
nent r\, the pole locations in the continued fraction expansion, 
are in excellent agreement in both cases. A smaller T] would 
have revealed the many (5-functions in the DMRG case ob- 
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FIG. 8: DMRG and exact results for G vs. V g in the case of 2 coupled 
non-interacting QD's. The DMRG results are obtained with L = 64 
and M = 300 (Al / =0.001). The exact results are for infinite leads. 
The cases (a) t" = 0.5 and (b) t" = 0.2 are investigated, with t' = 
0.4 in both cases. The results present the bonding and anti-bonding 
resonant tunneling peaks at ±t". 




FIG. 9: DMRG local density-of-states results for (a) 1QD and (b) 
2QD's compared to the exact (bulk) results. The parameters used 
are V s = U = 0, t' = 0.4 in both cases and t" = 0.5 in (b). In (a), a 
broadening imaginary component rj = 0.1 was introduced. In (b), 
r\ = 0.15. Smaller values of 77 would reveal the discrete nature of the 
LDOS obtained with DMRG on a finite L=64 system. 



tained using a finite chain with 64 sites. 



IV. ONE QUANTUM DOT 

In the previous sections, the method was introduced and 
tested for the case of noninteracting [/=0 electrons. But the 
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FIG. 10: J(t) / A V in the case of one quantum dot for different clus- 
ter lengths L. The parameters used are U = 1.0, if = 0.4, AV=0.001, 
and Vg = -0.5. As L increases, the conductance approaches the uni- 
tarity limit (2e 2 /h) due to the Kondo screening effect. 



main application of the technique is envisioned to occur in the 
presence of nontrivial Coulombic interactions (and eventually 
also adding phononic degrees of freedom). In this section, the 
case of a nonzero Hubbard coupling will be considered, focus- 
ing on the special case of one quantum dot. The Hamiltonian 
used was already discussed in previous sections. 



A. Results at intermediate values of U 

Figure ^3 contains our DMRG results for the current vs. 
time, for the case of t/ =1.0. Similar values of this coupling 
were extensively used in previous investigations^22iS2i^2i2iiS, 
and it is believed to lead to a Kondo cloud of a size amenable 
to numerical investigations (note that if U is very large, the ef- 
fective J between localized and mobile spins is reduced, and 
it is known that the cloud's size rapidly grows with decreas- 
ing J). The figure shows that the systematic behavior found 
in the noninteracting case survives the presence of a Coulomb 
interaction, namely the current develops plateaus that can be 
used to determine the conductance. For instance, this effect is 
clearly present for L=96 and 128, although for smaller sizes 
(shown for completeness) the maximum current is 10% to 
20% less than expected and one must be cautious with size 
effects. The value of the gate voltage is -U/2, which in the 
absence of the Kondo effect would locate the system in the 
conductance "valley" (implying a near zero conductance) be- 
tween the Coulomb blockade peaks at -U and 0. The figure 
shows that the method introduced in this paper is able to repro- 
duce the existence of a Kondo effect, since the conductance is 
actually very close to the ideal limit 2e 2 /h^^^, rather than 
being negligible. This is a highly nontrivial test that the pro- 
posed technique has passed. 

Results for other values of the gate voltage are shown in 
FigE] Moving away from V s =-U/2, the current is reduced, 
as discussed below in more detail. Note that the plateaus con- 
tain small oscillations as a function of time^ The procedure 
followed here to extract the current needed for the calculation 
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FIG. 11: J(t)/AV in the case of one interacting QD for different 
values of V g , and with Al/=0.001. The value of the conductance 
is obtained by averaging the current over an interval of time, corre- 
sponding to the steady state. The solid horizontal lines represent this 
time interval over which the average of the current is taken, and the 
value of the average. The parameters used are U = 1.0, t' = 0.4, L = 
128, and M = 350. 



of the conductance is to carry out averages over time, as shown 
in the figure, once the plateau region is reached. The size of 
the oscillations give an indication of the errors in the numeri- 
cal determination of the conductance, for a given lattice size. 
As L and M increases, the oscillations tend to disappear. 

Following the procedure sketched in Fig[^ the full con- 
ductance vs. V g was prepared for U equal to 1 and 2. The 
results are in Fig^] The shape of the curve is the expected 
one for the regime considered here: the intermediate values of 
U do not locate our investigation deep in the Kondo regime, 
with sharply defined integer charge at the dot, but more into 
the mixed-valence region. This can be deduced from the value 
of the dot charge vs. V g , also shown in Fig^l With increas- 
ing U, sharper charge steps are formed, but the Kondo cloud 
size increases, as discussed later. 



B. Results at large values of U 

There are cases where the technique gives results that are 
only qualitatively correct. While clearly further increases of 
the number of states and lattice sizes will improve the ac- 
curacy, it is important to judge if at least the essence of the 
physics has been captured by our proposed method. In Figll3l 
results for U=4 are shown. This is a representative of the 
"large" U regime, since it must be compared with t' (as op- 
posed to t=\) that is only 0.4 in this figure. Another indication 
that this U is large is in Fig lOf b) where a clear sharp quanti- 
zation of the charge inside the dot is observed. In this large- U 
regime, the DMRG conductance is shown in (a). Clearly, there 
is a substantial difference between the Friedel sum-rule esti- 
mation (which has the correct "box" shape in the gate voltage 
range -1 and 0) and the DMRG numbers. However, at least 
the fact that there must be a nonzero conductance at V g =-U/2 
was properly captured by the method. This example illus- 
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FIG. 12: Conductance G (AV=0.001) and the dot occupation (n d ) 
for one interacting QD. The circles show G obtained by averaging 
the current over an interval of time corresponding to the steady state, 
as shown in Fig ll II The squares show G obtained from (nd) using 
the Friedel sum rule (FSR). G has the shape of the expected Kondo or 
mixed- valence plateau centered at V s = —U/2. The feature would be 
sharper reducing t' . Results shown are: (a) U = 1.0, t' = 0.4, and (b) 
U = 2.0, t' = 0.5. In both cases L = 128 and M = 350. Note that the 
DMRG conductance results in (a) show a slightly better agreement 
with the FSR results. This is expected since the finite size effects are 
stronger for larger U. 



trates a case where size effects are important, due to the subtle 
rapid increase of the Kondo cloud with increasing U. How- 
ever, the qualitative results were captured, in spite of the fact 
that to reach a quantitative conclusion much larger sizes must 
be considered. This case is shown as a cautionary example to 
the readers, that must be alert of the limitations of the numer- 
ical methods. Note that while FSR results are excellent, for 
other arbitrary cases it would be unclear whether the Friedel 
sum-rule method is valid and, as a consequence, not always 
this procedure can be used. 



C. Improving the Convergence 

To carry out the investigations presented thus far in this and 
the previous section, the number of sites in the leads at left 
and right of the dot region have been chosen such that one lead 
has an even number and the other an odd number of sites (to 
refer to this case, the notation used below will be "odd-QDs- 
even"). While this should be irrelevant for very long chains, in 
practice this is important for the speed of convergence of the 
conductance calculations with increasing cluster size. For ex- 
ample, in the interacting case U^0, a Kondo or mixed-valence 
regime is expected where the spin at the dot couples with elec- 
trons at the Fermi level of the leads. This formation of the 
Kondo cloud occurs more efficiently on a finite-size lattice if 
already a zero energy level is available, as it occurs when one 
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FIG. 13: Conductance G (AV=0.001) and charge at the dot (n d ), 
for one interacting QD in the case of large U. The circles show 
G obtained from DMRG current procedure outlines in this paper, 
while the squares show G obtained from (rid) using the Friedel sum- 
rule. The finite-size effects are obvious here, since the results are 
halfway between the expected Kondo plateau (properly reproduced 
by the FSR procedure) and the Coulomb blockade peaks. This case 
is shown as an illustration of important size effects in some limits. 
The parameters used are U = 4.0, t' = 0.4, L = 128, and M = 300. 



of the leads has an odd number of sites. That this improves 
the rate of convergence with increasing lattice size is shown 
in Fig^J where the odd-lQD-even case is contrasted to the 
even-lQD-even case, where in both leads the number of sites 
is even. Clearly, the odd-lQD-even case approaches the ideal 
limit 2e 2 lh faster than using leads with an even number of 
sites, and it is recommended to be used in future investiga- 
tions. The figure also shows that eventually with sufficiently 
large systems, both combinations will reach the same ideal 
result, as expected. The remaining possibility (odd-lQD-odd) 
was also investigated (not shown). For reasons that remain 
to be analyzed further, the convergence in this case is not as 
good as in the odd- lQD-even case. As a consequence, empiri- 
cally it is clear that the combination even-QDs-odd is the most 
optimal to speed up the size convergence of the calculation. 

Another method to improve the size convergence was 
tested. Following Ref. [251 the finite-size effects can be re- 
duced by using "damped boundary conditions" (DBC). The 
hoppings in the Md links at the boundaries are reduced using 
the formulas —td : —td 2 , —td M °, where d<l. Mr> has to 
be chosen such that the damping occurs far enough from the 
central region. Figure IT31 shows the finite size scaling using 
odd-QD-even clusters and the same parameters as in Fig ll4l 
both with DBC and the standard open boundary conditions 
(OBC) used in the rest of the paper. The latter indeed im- 
proves the convergence. While in the present study, OBC 
were used in most of the manuscript to keep the simplicity 
in the presentation and reduce the number of parameters, the 
use of DBC is recommended for cases where size effects are 
strong. 
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FIG. 14: Finite-size scaling of the conductance G at V g = -U /2 for 
the odd-lQD-even (circles) and even-lQD-even clusters (triangles) 
setups. Note that in both cases G converges to 1 in the bulk limit. 
However, the odd-lQD-even cluster converges faster, which makes 
it the most useful for practical calculations. The parameters used are 
U = 1.0 and £'=0.4. 
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FIG. 15: Finite size scaling of G for odd-QD-even clusters using 
OBC (circles) and DBC (triangles), d — 0.5 was used. Note that 
L — 64 cluster with DBC gives better results than a L=128 cluster 
with OBC. The same parameters are used as in Fie ll4l 



D. Influence of Magnetic Fields 

To fully confirm that our investigations in the Kondo/mixed 
valence regime have captured the essence of the problem, 
namely the formation of a Kondo cloud with antiferromag- 
netic coupling between the spins at the leads and the dot, 
investigations including magnetic fields are necessary. In 
FigtHD it is shown how the the Kondo plateau in the con- 
ductance evolves with increasing magnetic field. As ex- 
pected from previous investigations, including results ob- 
tained using very different techniques such as the Lanczos 
method followed by a Dyson-equation embedding procedure 
(ED+DE) 27 , the conductance broad peak splits with increas- 
ing magnetic field B. At large B, two peaks are observed 
at -U and 0, as it occurs also in the high temperature regime 
where only Coulomb blockade effects are present. 
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FIG. 16: Conductance G vs. gate voltage Vg in the case of one inter- 
acting QD (U = 2.0, t' = 0.5) for different values of the magnetic field 
B (and AV=0.001). For B = 0, a Kondo plateau is obtained, cen- 
tered at V g = —U /2. As B increases, the Kondo effect is suppressed, 
and for moderate B, two Coulomb blockade peaks are observed at 
V g = — 17 and V s = 0, as expected. 
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FIG. 17: J(t) /A V for two coupled QD's at V g = -U/2 for different 
values of t" /t' 2 (and AV=0.001). As in the case of one interacting 
QD, the conductance is obtained by averaging the steady-state cur- 
rent over the indicated intervals. The parameters used are U =1.0, t' 
= 0.5, L = 127, and M = 350. 



V. TWO COUPLED QUANTUM DOTS 

The method discussed in this paper is general, and in prin- 
ciple it can be implemented for a variety of complicated ge- 
ometries and couplings in the interactive region between the 
leads. Thus, it is important to confirm that the method will 
keep its validity going beyond the one quantum dot case. In 
this section, the case of two dots will be studied. Systems with 
two quantum dots are believed to be understood theoretically 
and, as a consequence, our numerical data can be contrasted 
against robust results in the literature. Cases involving more 
dots2& are still not fully understood, and their analysis will be 
postponed for future investigations. In Fig^] the current vs. 
time is shown for two dots. The Hamiltonian for this case was 
already defined in previous sections. For a fixed t' , increas- 
ing the amplitude of the direct hopping between the dots t" 
amounts to isolating the two dots system from the rest. As 
a consequence, the current is expected to decrease, and the 
method indeed reproduces this effect, as shown in the figure. 
The same physics is obtained reducing t', at fixed t". In fact, 
previous studies 35 have shown that the conductance only de- 
pends on t"/t' 2 , and this has been verified using our method. 

The conductance of a system with two dots in series will de- 
crease with increasing t" It' 2 (at large t" It 12 ') due to the decou- 
pling of the two-dots system into a small two-sites molecule, 
as already discussed. But this conductance will also be very 
small at small t" when the tunneling from one dot to the next 
is nearly cutoff. As a consequence, the conductance vs. t" It' 2 
is known to present a peak at intermediate values. In Figll 81 
the slave-boson mean-field technique (SBMFT) predictions 
for this case obtained in previous investigations 35 are shown 
together with our results. The agreement is fairly reasonable, 
providing further support that the method discussed here can 
handle systems where there are competing tendencies, beyond 
the one quantum-dot case. 
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FIG. 18: Conductance G as a function of t"/t' 2 , at V g = -(7/2 
(A V=0.001). In this regime, G is determined by the competition be- 
tween the Kondo correlation of each dot with the neighboring leads 
and the antiferromagnetic correlation between the two dots. The cir- 
cles represent our DMRG results obtained with L = 127 and M = 
350. The solid line is the plot of the functional form obtained by 
Georges and Meir using SBMFT^ 



VI. CONCLUSIONS 

In this manuscript, a method was proposed and tested to 
calculate the conductance of small (nanoscale) strongly corre- 
lated systems modeled by tight-binding Hamiltonians. The 
approach is based on the adaptive time-dependent DMRG 
method, and it was shown to work properly for non-interacting 
systems and also in the cases of one and two quantum dots. 
Besides the finite size effects, discussed in the text as well, 
there are no other severe limitations to handle complex in- 
teracting models with arbitrary couplings. The method is a 
complement to DFT calculations in the nanoscopic context. 
Further improvements of the technique must consider temper- 
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ature and inelastic effects. 

It is concluded that the semi-quantitative analysis of trans- 
port in models of strongly correlated nanosystems appears un- 
der reach, and much progress is expected from the application 
of the technique presented here to realistic Hamiltonians for 
small molecules and arrays of quantum dots. The remarkable 
cross-fertilization between modeling, simulation, and experi- 
ments existing in bulk strongly correlated materials, such as 
transition-metal oxides, can be repeated in a variety of inter- 
esting systems at the nanoscale. 
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